## R Script Output -------------------------------------------------------------
# Figure 3: Effect Heterogeneity: Nativism Mechanisms.
# Figure 4: Effect Heterogeneity: Pocketbook Voting Mechanism
# Appendix Figure D.2: Heterogeneous Treatment Effects: Joint Model Results
# Appendix Table D.4: Heterogeneous Treatment Effects: Nativism Hypotheses
# Appendix Table D.5: Heterogeneous Treatment Effects: Pocketbook Voting Hypotheses


## Instructions ----------------------------------------------------------------
# Step 1: Adjust MAIN_DIR to where README.txt is located
# Step 2: Run entire script


## setup -----------------------------------------------------------------------
# clean slate
rm(list = ls())
date()

# load packages
pkg <- c("tidyverse",
         "RColorBrewer", 
         "gridExtra",
         "tools", 
         "lfe", 
         "stargazer",
         "multcomp",
         "pBrackets")

lapply(pkg, require, character.only = TRUE)

# set main directory
MAIN_DIR <- "~/Dropbox/Research/ISQ-frei-replication/"


## load data -------------------------------------------------------------------
load(file = paste(MAIN_DIR, "data-merge-county.RData", sep = ""))

# set data frames
df <- county.df

df.treat.1 <- df %>% 
  filter(year == 2012 | year == 2016)

df.treat.2 <- df %>% 
  filter(year == 2012 | year == 2020)

df.placebo <- df %>% 
  filter(year == 2008 | year == 2012)


## set variable and variable name correspondence -------------------------------
var.df <- tibble(var = c("anticorrupt",
                         "placebo_anticorrupt_2012",
                         
                         "anticorrupt:cn_ug_sqmi_cat_2012_bi_high",
                         "anticorrupt:cn_ug_sqmi_cat_2012_tri_medium",
                         "anticorrupt:cn_ug_sqmi_cat_2012_tri_high",
                         
                         "anticorrupt:cn_g_sqmi_cat_2012_biCN-G-high",
                         "anticorrupt:cn_g_sqmi_cat_2012_triCN-G-medium",
                         "anticorrupt:cn_g_sqmi_cat_2012_triCN-G-high",
                         
                         "anticorrupt:ind_ug_sqmi_cat_2012_biIND-UG-high",
                         "anticorrupt:ind_ug_sqmi_cat_2012_triIND-UG-medium",
                         "anticorrupt:ind_ug_sqmi_cat_2012_triIND-UG-high",
                         
                         "placebo_anticorrupt_2012:cn_ug_sqmi_cat_2012_bi_high",
                         "placebo_anticorrupt_2012:cn_ug_sqmi_cat_2012_tri_medium",
                         "placebo_anticorrupt_2012:cn_ug_sqmi_cat_2012_tri_high",
                         
                         "placebo_anticorrupt_2012:cn_g_sqmi_cat_2012_biCN-G-high",
                         "placebo_anticorrupt_2012:cn_g_sqmi_cat_2012_triCN-G-medium",
                         "placebo_anticorrupt_2012:cn_g_sqmi_cat_2012_triCN-G-high",
                         
                         "placebo_anticorrupt_2012:ind_ug_sqmi_cat_2012_biIND-UG-high",
                         "placebo_anticorrupt_2012:ind_ug_sqmi_cat_2012_triIND-UG-medium",
                         "placebo_anticorrupt_2012:ind_ug_sqmi_cat_2012_triIND-UG-high",
                         
                         # controls
                         "share_pop_county_female",
                         "share_pop_5_17_county",
                         "share_pop_18_24_county",
                         "share_pop_25_34_county",
                         "share_pop_35_44_county",
                         "share_pop_45_54_county",
                         "share_pop_55_64_county",
                         "share_pop_65_county",
                         
                         "share_pop_white_county",
                         "share_pop_black_county",
                         "share_pop_hispanic_county",
                         "share_pop_asian_county",
                         "share_pop_cn_county",
                         
                         "share_foreign_born_county",
                         "share_edu_ba_above_county",
                         "share_enroll_county",
                         "ln_pop_density_county",
                         "effective_tax_county", 
                         "ipw_county",
                         "employ_rate_county", 
                         "median_household_income_county",
                         "vacancy_county",
                         
                         # trending effects of controls
                         "anticorrupt:share_pop_county_female",
                         "anticorrupt:share_pop_5_17_county",
                         "anticorrupt:share_pop_18_24_county",
                         "anticorrupt:share_pop_25_34_county",
                         "anticorrupt:share_pop_35_44_county",
                         "anticorrupt:share_pop_45_54_county",
                         "anticorrupt:share_pop_55_64_county",
                         "anticorrupt:share_pop_65_county",
                         
                         "anticorrupt:share_pop_white_county",
                         "anticorrupt:share_white_2012_cat_medium",
                         "anticorrupt:share_white_2012_cat_high",
                         "anticorrupt:share_pop_black_county",
                         "anticorrupt:share_pop_hispanic_county",
                         "anticorrupt:share_pop_asian_county",
                         "anticorrupt:share_pop_cn_county",
                         
                         "anticorrupt:share_foreign_born_county",
                         "anticorrupt:share_edu_ba_above_county",
                         "anticorrupt:share_edu_ba_above_county_2012_tri_medium",
                         "anticorrupt:share_edu_ba_above_county_2012_tri_high",
                         "anticorrupt:share_enroll_county",
                         "anticorrupt:ln_pop_density_county",
                         "anticorrupt:effective_tax_county", 
                         "anticorrupt:ipw_county",
                         "anticorrupt:employ_rate_county", 
                         "anticorrupt:median_household_income_county",
                         "anticorrupt:vacancy_county",
                         "anticorrupt:vacancy_county_2012_tri_medium",
                         "anticorrupt:vacancy_county_2012_tri_high",
                         
                         "anticorrupt:share_homeowner_county_2012_cat_medium",
                         "anticorrupt:share_homeowner_county_2012_cat_high",
                         
                         # trending effects of controls in placebo period
                         "placebo_anticorrupt_2012:share_pop_county_female",
                         "placebo_anticorrupt_2012:share_pop_5_17_county",
                         "placebo_anticorrupt_2012:share_pop_18_24_county",
                         "placebo_anticorrupt_2012:share_pop_25_34_county",
                         "placebo_anticorrupt_2012:share_pop_35_44_county",
                         "placebo_anticorrupt_2012:share_pop_45_54_county",
                         "placebo_anticorrupt_2012:share_pop_55_64_county",
                         "placebo_anticorrupt_2012:share_pop_65_county",
                         
                         "placebo_anticorrupt_2012:share_pop_white_county",
                         "placebo_anticorrupt_2012:share_pop_black_county",
                         "placebo_anticorrupt_2012:share_pop_hispanic_county",
                         "placebo_anticorrupt_2012:share_pop_asian_county",
                         "placebo_anticorrupt_2012:share_pop_cn_county",
                         
                         "placebo_anticorrupt_2012:share_foreign_born_county",
                         "placebo_anticorrupt_2012:share_edu_ba_above_county",
                         "placebo_anticorrupt_2012:share_enroll_county",
                         "placebo_anticorrupt_2012:ln_pop_density_county",
                         "placebo_anticorrupt_2012:effective_tax_county", 
                         "placebo_anticorrupt_2012:ipw_county",
                         "placebo_anticorrupt_2012:employ_rate_county", 
                         "placebo_anticorrupt_2012:median_household_income_county",
                         "placebo_anticorrupt_2012:vacancy_county",
                         
                         # triple interactions
                         "anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium",
                         "anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high",
                         "anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium",
                         "anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high",
                         "anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium",
                         "anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high",
                         "anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium",
                         "anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high"
),
var_name = c("2013 Anti-corruption Campaign (Dummy)",
             "Placebo Anti-corruption (2012)",
             
             "2013 Anti-corruption X CN UG in 2012 (Dummy, High)",
             "2013 Anti-corruption X CN UG in 2012 (Tri., Medium)",
             "2013 Anti-corruption X CN UG in 2012 (Tri., High)",
             
             "2013 Anti-corruption X CN G in 2012 (Dummy, High)",
             "2013 Anti-corruption X CN G in 2012 (Tri., Medium)",
             "2013 Anti-corruption X CN G in 2012 (Tri., High)",
             
             "2013 Anti-corruption X IND UG in 2012 (Dummy, High)",
             "2013 Anti-corruption X IND UG in 2012 (Tri., Medium)",
             "2013 Anti-corruption X IND UG in 2012 (Tri., High)",
             
             "Placebo Anti-corruption (2012) X CN UG in 2012 (Dummy, High)",
             "Placebo Anti-corruption (2012) X CN UG in 2012 (Tri., Medium)",
             "Placebo Anti-corruption (2012) X CN UG in 2012 (Tri., High)",
             
             "Placebo Anti-corruption (2012) X CN G in 2012 (Dummy, High)",
             "Placebo Anti-corruption (2012) X CN G in 2012 (Tri., Medium)",
             "Placebo Anti-corruption (2012) X CN G in 2012 (Tri., High)",
             
             "Placebo Anti-corruption (2012) X IND UG in 2012 (Dummy, High)",
             "Placebo Anti-corruption (2012) X IND UG in 2012 (Tri., Medium)",
             "Placebo Anti-corruption (2012) X IND UG in 2012 (Tri., High)",
             
             # controls
             "Share of Female Pop.",
             "Share of Pop. 5--17",
             "Share of Pop. 18--24",
             "Share of Pop. 25--34",
             "Share of Pop. 35--44",
             "Share of Pop. 45--54",
             "Share of Pop. 55--64",
             "Share of Pop. 65--",
             
             "Share of White Pop.",
             "Share of Black Pop.",
             "Share of Hispanic Pop.",
             "Share of Asian Pop.",
             "Share of Chinese Pop.",
             
             "Share of Foreign-Born Pop.",
             "Share of Pop. with BA Degree or Above",
             "Share of Pop. Enrolled in College or Above",
             "Population Density (log)",
             "Effective Tax Rate (\\%)",
             "Trade Exposure (IPW)",
             "Employment Rate",
             "Median Household Income (10,000)",
             "Share of Vacant Houses",
             
             # trending effects of controls
             "2013 Anti-corruption X Share of Female Pop.",
             "2013 Anti-corruption X Share of Pop. 5--17",
             "2013 Anti-corruption X Share of Pop. 18--24",
             "2013 Anti-corruption X Share of Pop. 25--34",
             "2013 Anti-corruption X Share of Pop. 35--44",
             "2013 Anti-corruption X Share of Pop. 45--54",
             "2013 Anti-corruption X Share of Pop. 55--64",
             "2013 Anti-corruption X Share of Pop. 65--",
             
             "2013 Anti-corruption X Share of White Pop.",
             "2013 Anti-corruption X 2012 White Pop. (Medium)",
             "2013 Anti-corruption X 2012 White Pop. (Large)",
             "2013 Anti-corruption X Share of Black Pop.",
             "2013 Anti-corruption X Share of Hispanic Pop.",
             "2013 Anti-corruption X Share of Asian Pop.",
             "2013 Anti-corruption X Share of Chinese Pop.",
             
             "2013 Anti-corruption X Share of Foreign-Born Pop.",
             "2013 Anti-corruption X Share of Pop. with BA Degree or Above",
             "2013 Anti-corruption X 2012 College Educated Pop. (Medium)",
             "2013 Anti-corruption X 2012 College Educated Pop. (Large)",
             "2013 Anti-corruption X Share of Pop. Enrolled in College or Above",
             "2013 Anti-corruption X Population Density (log)",
             "2013 Anti-corruption X Effective Tax Rate (\\%)",
             "2013 Anti-corruption X Trade Exposure (IPW)",
             "2013 Anti-corruption X Employment Rate",
             "2013 Anti-corruption X Median Household Income (10,000)",
             "2013 Anti-corruption X Share of Vacant Houses",
             "2013 Anti-corruption X 2012 Vacancy Rate (Medium)",
             "2013 Anti-corruption X 2012 Vacancy Rate (High)",
             
             "2013 Anti-corruption X 2012 Homeownership Rate (Medium)", 
             "2013 Anti-corruption X 2012 Homeownership Rate (High)", 
             
             # trending effects of controls in placebo period
             "Placebo Anti-corruption (2012) X Share of Female Pop.",
             "Placebo Anti-corruption (2012) X Share of Pop. 5--17",
             "Placebo Anti-corruption (2012) X Share of Pop. 18--24",
             "Placebo Anti-corruption (2012) X Share of Pop. 25--34",
             "Placebo Anti-corruption (2012) X Share of Pop. 35--44",
             "Placebo Anti-corruption (2012) X Share of Pop. 45--54",
             "Placebo Anti-corruption (2012) X Share of Pop. 55--64",
             "Placebo Anti-corruption (2012) X Share of Pop. 65--",
             
             "Placebo Anti-corruption (2012) X Share of White Pop.",
             "Placebo Anti-corruption (2012) X Share of Black Pop.",
             "Placebo Anti-corruption (2012) X Share of Hispanic Pop.",
             "Placebo Anti-corruption (2012) X Share of Asian Pop.",
             "Placebo Anti-corruption (2012) X Share of CN Pop.",
             
             "Placebo Anti-corruption (2012) X Share of Foreign-Born Pop.",
             "Placebo Anti-corruption (2012) X Share of Pop. with BA Degree or Above",
             "Placebo Anti-corruption (2012) X Share of Pop. Enrolled in College or Above",
             "Placebo Anti-corruption (2012) X Population Density (log)",
             "Placebo Anti-corruption (2012) X Effective Tax Rate (\\%)",
             "Placebo Anti-corruption (2012) X Trade Exposure (IPW)",
             "Placebo Anti-corruption (2012) X Employment Rate",
             "Placebo Anti-corruption (2012) X Median Household Income (10,000)",
             "Placebo Anti-corruption (2012) X Share of Vacant Houses",
             
             # triple interactions
             "2013 Anti-corruption X CN UG in 2012 (Dummy, High) X 2012 White Population (Medium)",
             "2013 Anti-corruption X CN UG in 2012 (Dummy, High) X 2012 White Population (Large)",
             "2013 Anti-corruption X CN UG in 2012 (Dummy, High) X 2012 College Educated Pop. (Medium)",
             "2013 Anti-corruption X CN UG in 2012 (Dummy, High) X 2012 College Educated Pop. (Large)",
             "2013 Anti-corruption X CN UG in 2012 (Dummy, High) X 2012 Homeownership Rate (Medium)",
             "2013 Anti-corruption X CN UG in 2012 (Dummy, High) X 2012 Homeownership Rate (High)",
             "2013 Anti-corruption X CN UG in 2012 (Dummy, High) X 2012 Vacancy Rate (Medium)",
             "2013 Anti-corruption X CN UG in 2012 (Dummy, High) X 2012 Vacancy Rate (High)"
))


## function to replace variable names ------------------------------------------
replaceVarName <- function(var.vec, var.df){
  # Prepare output vector
  out.vec <- rep(NA, length(var.vec))
  matches <- match(var.vec, var.df$var)
  out.vec <- var.df[matches,]$var_name
  
  if(any(is.na(out.vec))){
    warning(paste("Variable concordence missing: ", 
                  paste(var.vec[is.na(out.vec)], collapse = ", "), 
                  sep = ""))
  } else{
    #print("All variables successfully converted")
  }
  
  return(out.vec)
}


## Table D.4, Column 1 ---------------------------------------------------------
ddd.12.16.w.f <- felm(dem_share ~ 
                        anticorrupt + 
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high + 
                        
                        anticorrupt:share_white_2012_cat_medium + 
                        anticorrupt:share_white_2012_cat_high + 
                        
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high +
                        
                        # controls
                        cn_g_sqmi_cat_2012_bi:anticorrupt +
                        ind_ug_sqmi_cat_2012_bi:anticorrupt +
                        
                        # demographics
                        share_pop_county_female * anticorrupt +
                        
                        share_pop_5_17_county * anticorrupt +
                        share_pop_18_24_county * anticorrupt +
                        share_pop_25_34_county * anticorrupt +
                        share_pop_35_44_county * anticorrupt +
                        share_pop_45_54_county * anticorrupt +
                        share_pop_55_64_county * anticorrupt +
                        share_pop_65_county * anticorrupt + 
                        
                        share_pop_white_county * anticorrupt + 
                        share_pop_black_county * anticorrupt + 
                        share_pop_hispanic_county * anticorrupt + 
                        share_pop_asian_county * anticorrupt + 
                        share_pop_cn_county * anticorrupt + 
                        
                        share_foreign_born_county * anticorrupt + 
                        
                        share_edu_ba_above_county * anticorrupt +
                        share_enroll_county * anticorrupt + 
                        ln_pop_density_county * anticorrupt + 
                        
                        effective_tax_county * anticorrupt + 
                        
                        ipw_county * anticorrupt + 
                        
                        employ_rate_county * anticorrupt + 
                        median_household_income_county * anticorrupt + 
                        vacancy_county * anticorrupt 
                      | county_fips | 0 | county_fips, 
                      keepX = TRUE,
                      data = df.treat.1)
summary(ddd.12.16.w.f)
n_distinct(ddd.12.16.w.f$fe$county_fips)


## Table D.4, Column 2 ---------------------------------------------------------
ddd.12.20.w.f <- felm(dem_share ~ 
                        anticorrupt + 
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high + 
                        
                        anticorrupt:share_white_2012_cat_medium + 
                        anticorrupt:share_white_2012_cat_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high +
                        
                        # controls
                        cn_g_sqmi_cat_2012_bi:anticorrupt +
                        ind_ug_sqmi_cat_2012_bi:anticorrupt +
                        
                        # demographics
                        share_pop_county_female * anticorrupt +
                        
                        share_pop_5_17_county * anticorrupt +
                        share_pop_18_24_county * anticorrupt + 
                        share_pop_25_34_county * anticorrupt +
                        share_pop_35_44_county * anticorrupt +
                        share_pop_45_54_county * anticorrupt +
                        share_pop_55_64_county * anticorrupt +
                        share_pop_65_county * anticorrupt + 
                        
                        share_pop_white_county * anticorrupt + 
                        share_pop_black_county * anticorrupt + 
                        share_pop_hispanic_county * anticorrupt + 
                        share_pop_asian_county * anticorrupt + 
                        share_pop_cn_county * anticorrupt + 
                        
                        share_foreign_born_county * anticorrupt + 
                        
                        share_edu_ba_above_county * anticorrupt +
                        share_enroll_county * anticorrupt + 
                        ln_pop_density_county * anticorrupt + 
                        
                        effective_tax_county * anticorrupt + 
                        
                        ipw_county * anticorrupt + 
                        
                        employ_rate_county * anticorrupt + 
                        median_household_income_county * anticorrupt + 
                        vacancy_county * anticorrupt 
                      | county_fips | 0 | county_fips, 
                      keepX = TRUE,
                      data = df.treat.2)
summary(ddd.12.20.w.f)
n_distinct(ddd.12.20.w.f$fe$county_fips)


## Table D.4, Column 3 ---------------------------------------------------------
ddd.12.16.e.f <- felm(dem_share ~ 
                        anticorrupt + 
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high + 
                        
                        anticorrupt:share_edu_ba_above_county_2012_tri_medium + 
                        anticorrupt:share_edu_ba_above_county_2012_tri_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high +
                        
                        # controls
                        cn_g_sqmi_cat_2012_bi:anticorrupt +
                        ind_ug_sqmi_cat_2012_bi:anticorrupt +
                        
                        # demographics
                        share_pop_county_female * anticorrupt +
                        
                        share_pop_5_17_county * anticorrupt +
                        share_pop_18_24_county * anticorrupt + 
                        share_pop_25_34_county * anticorrupt +
                        share_pop_35_44_county * anticorrupt +
                        share_pop_45_54_county * anticorrupt +
                        share_pop_55_64_county * anticorrupt +
                        share_pop_65_county * anticorrupt + 
                        
                        share_pop_white_county * anticorrupt + 
                        share_pop_black_county * anticorrupt + 
                        share_pop_hispanic_county * anticorrupt + 
                        share_pop_asian_county * anticorrupt + 
                        share_pop_cn_county * anticorrupt +
                        
                        share_foreign_born_county * anticorrupt + 
                        
                        share_edu_ba_above_county * anticorrupt +
                        share_enroll_county * anticorrupt + 
                        ln_pop_density_county * anticorrupt + 
                        
                        effective_tax_county * anticorrupt + 
                        
                        ipw_county * anticorrupt + 
                        
                        employ_rate_county * anticorrupt + 
                        median_household_income_county * anticorrupt + 
                        vacancy_county * anticorrupt
                      | county_fips | 0 | county_fips, 
                      keepX = TRUE,
                      data = df.treat.1)
summary(ddd.12.16.e.f)
n_distinct(ddd.12.16.e.f$fe$county_fips)


## Table D.4, Column 4 ---------------------------------------------------------
ddd.12.20.e.f <- felm(dem_share ~ 
                        anticorrupt + 
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high + 
                        
                        anticorrupt:share_edu_ba_above_county_2012_tri_medium + 
                        anticorrupt:share_edu_ba_above_county_2012_tri_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high +
                        
                        # controls
                        cn_g_sqmi_cat_2012_bi:anticorrupt +
                        ind_ug_sqmi_cat_2012_bi:anticorrupt +
                        
                        # demographics
                        share_pop_county_female * anticorrupt +
                        
                        share_pop_5_17_county * anticorrupt +
                        share_pop_18_24_county * anticorrupt + 
                        share_pop_25_34_county * anticorrupt +
                        share_pop_35_44_county * anticorrupt +
                        share_pop_45_54_county * anticorrupt +
                        share_pop_55_64_county * anticorrupt +
                        share_pop_65_county * anticorrupt + 
                        
                        share_pop_white_county * anticorrupt + 
                        share_pop_black_county * anticorrupt + 
                        share_pop_hispanic_county * anticorrupt + 
                        share_pop_asian_county * anticorrupt + 
                        share_pop_cn_county * anticorrupt + 
                        
                        share_foreign_born_county * anticorrupt + 
                        
                        share_edu_ba_above_county * anticorrupt +
                        share_enroll_county * anticorrupt + 
                        ln_pop_density_county * anticorrupt + 
                        
                        effective_tax_county * anticorrupt + 
                        
                        ipw_county * anticorrupt + 
                        
                        employ_rate_county * anticorrupt + 
                        median_household_income_county * anticorrupt + 
                        vacancy_county * anticorrupt 
                      | county_fips | 0 | county_fips, 
                      keepX = TRUE,
                      data = df.treat.2)
summary(ddd.12.20.e.f)
n_distinct(ddd.12.20.e.f$fe$county_fips)


## Table D.4: Heterogeneous Treatment Effects: Nativism Hypotheses -------------
# set list of models
did.nativism <- list(ddd.12.16.w.f,
                     ddd.12.20.w.f,
                     ddd.12.16.e.f,
                     ddd.12.20.e.f)

# set list of clustered SEs
did.nativism.cse <- list(ddd.12.16.w.f$cse,
                         ddd.12.20.w.f$cse,
                         ddd.12.16.e.f$cse,
                         ddd.12.20.e.f$cse)
  

# set model names
mod.names.nativism <- c( "12-16-full",
                         "12-20-full",
                         "12-16-full",
                         "12-20-full")

# set variable order
var.order.nativism <- c("^anticorrupt$",
                        "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high$",
                        
                        "^anticorrupt:share_white_2012_cat_medium$",
                        "^anticorrupt:share_white_2012_cat_high$",
                        "^anticorrupt:share_edu_ba_above_county_2012_tri_medium$",
                        "^anticorrupt:share_edu_ba_above_county_2012_tri_high$",
                        
                        "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium$",
                        "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high$",
                        "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium$",
                        "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high$",
                        
                        # controls
                        "^anticorrupt:cn_g_sqmi_cat_2012_biCN-G-high$",
                        "^anticorrupt:ind_ug_sqmi_cat_2012_biIND-UG-high$",
                        
                        "^share_pop_county_female$",
                        "^share_pop_5_17_county$",
                        "^share_pop_18_24_county$",
                        "^share_pop_25_34_county$",
                        "^share_pop_35_44_county$",
                        "^share_pop_45_54_county$",
                        "^share_pop_55_64_county$",
                        "^share_pop_65_county$",
                        
                        "^share_pop_white_county$",
                        "^share_pop_black_county$",
                        "^share_pop_hispanic_county$",
                        "^share_pop_asian_county$",
                        "^share_pop_cn_county$",
                        
                        "^share_foreign_born_county$",
                        "^share_edu_ba_above_county$",
                        "^share_enroll_county$",
                        "^ln_pop_density_county$",
                        "^effective_tax_county$", 
                        "^ipw_county$",
                        "^employ_rate_county$", 
                        "^median_household_income_county$",
                        "^vacancy_county$",
                        
                        # trending effects of controls
                        "^anticorrupt:share_pop_county_female$",
                        "^anticorrupt:share_pop_5_17_county$",
                        "^anticorrupt:share_pop_18_24_county$",
                        "^anticorrupt:share_pop_25_34_county$",
                        "^anticorrupt:share_pop_35_44_county$",
                        "^anticorrupt:share_pop_45_54_county$",
                        "^anticorrupt:share_pop_55_64_county$",
                        "^anticorrupt:share_pop_65_county$",
                        
                        "^anticorrupt:share_pop_white_county$",
                        "^anticorrupt:share_pop_black_county$",
                        "^anticorrupt:share_pop_hispanic_county$",
                        "^anticorrupt:share_pop_asian_county$",
                        "^anticorrupt:share_pop_cn_county$",
                        
                        "^anticorrupt:share_foreign_born_county$",
                        "^anticorrupt:share_edu_ba_above_county$",
                        "^anticorrupt:share_enroll_county$",
                        "^anticorrupt:ln_pop_density_county$",
                        "^anticorrupt:effective_tax_county$", 
                        "^anticorrupt:ipw_county$",
                        "^anticorrupt:employ_rate_county$", 
                        "^anticorrupt:median_household_income_county$",
                        "^anticorrupt:vacancy_county$")

# set var label
var.label.nativism <- str_replace_all(var.order.nativism, "\\^", "")
var.label.nativism <- str_replace_all(var.label.nativism, "\\$", "")

# save
sink(file.path(MAIN_DIR, "Table-D4.txt"))
#sink(file.path(MAIN_DIR, "Table-D4.tex"))
stargazer(did.nativism,
          #type = "latex",
          type = "text",
          title = "Heterogeneous Treatment Effects: Nativism Hypotheses",
          dep.var.labels = "Two-Party Democratic Vote Share",
          dep.var.labels.include = TRUE,
          column.labels = mod.names.nativism,
          order = var.order.nativism,
          covariate.labels = replaceVarName(var.vec = var.label.nativism,
                                            var.df = var.df),
          add.lines = list(c("Fixed Effects: County",
                             rep("Yes", 4)),
                           c("Number of Counties",
                             formatC(n_distinct(ddd.12.16.w.f$fe$county_fips), big.mark = ","),
                             formatC(n_distinct(ddd.12.20.w.f$fe$county_fips), big.mark = ","),
                             formatC(n_distinct(ddd.12.16.e.f$fe$county_fips), big.mark = ","),
                             formatC(n_distinct(ddd.12.20.e.f$fe$county_fips), big.mark = ",")
                           )),
          label = "tb:did-pres-moderation-nativism",
          single.row = FALSE,
          se = did.nativism.cse,
          df = FALSE,
          font.size = "tiny",
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          star.char = c("+", "*", "**", "***"),
          notes.append = FALSE)
sink()


## Figure 3: Effect Heterogeneity: Nativism Mechanisms -------------------------
# extract DiD results in list
fm.list.nativism <- list("Initial White Population\n" = summary(glht(ddd.12.16.w.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial White Population\n" = summary(glht(ddd.12.16.w.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium = 0"))),
                         "Initial White Population\n" = summary(glht(ddd.12.16.w.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high = 0"))),
                         
                         "Initial White Population\n" = summary(glht(ddd.12.20.w.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial White Population\n" = summary(glht(ddd.12.20.w.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium = 0"))),
                         "Initial White Population\n" = summary(glht(ddd.12.20.w.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high = 0"))),
                         
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.16.e.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.16.e.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium = 0"))),
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.16.e.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high = 0"))),
                         
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.20.e.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.20.e.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium = 0"))),
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.20.e.f, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high = 0")))
)

# extract key coefs and clustered SEs
coef.df.nativism <- map_df(fm.list.nativism, function(x){ 
  out <- tibble(var = names(coef(x)),
                coef = coef(x),
                cse = sqrt(diag(vcov(x))))
}, .id = "model") 

# add vars
var.name.nativism <- c("Small",
                       "Medium",
                       "Large",
                       "Small",
                       "Medium",
                       "Large")

var.period.nativism <- c(rep("2012 vs. 2016", 3),
                         rep("2012 vs. 2020", 3))

# combine vars
coef.df.nativism <- cbind(coef.df.nativism, 
                          var.name.nativism,
                          var.period.nativism)

# set alpha level
alpha <- 0.05

# compute C.I.
coef.df.nativism <- coef.df.nativism %>%
  mutate(model = factor(model, 
                        levels = c("Initial White Population\n",
                                   "Initial College-Educated Population\n")),
         multiplier = qnorm(1 - alpha / 2),
         coef = as.numeric(coef),
         cse = as.numeric(cse),
         lwr = coef - multiplier * cse,
         upr = coef + multiplier * cse,
         mod_var = paste(model, var.name.nativism, var.period.nativism, sep = "-"))

# set var order
coef.df.nativism <- coef.df.nativism %>%
  mutate(mod_var = factor(mod_var, 
                          levels = coef.df.nativism$mod_var))

# set text df for annotations
txt.shift.nativism <- 0.002
txt.pos.nativism <- tibble(model = coef.df.nativism$model,
                           var = as.character(coef.df.nativism$var.name),
                           mod_var = coef.df.nativism$mod_var,
                           coef = coef.df.nativism$lwr - txt.shift.nativism,
                           lwr = coef.df.nativism$lwr,
                           upr = coef.df.nativism$upr)

# set parameters
axis.title.size <- 16

# add curly brackets
bracketsGrob <- function(...){
  l <- list(...)
  e <- new.env()
  e$l <- l
  grid:::recordGrob(  {
    do.call(grid.brackets, l)
  }, e)
}

brack.y <- 0.87
b1 <- bracketsGrob(0.05, brack.y, 0.45, brack.y, 
                   h = 0.05, lwd = 1,
                   col = "black")
b2 <- bracketsGrob(0.55, brack.y, 0.95, brack.y, 
                   h = 0.05, lwd = 1,
                   col = "black")

# plot
plot.coef.nativism <- ggplot(data = coef.df.nativism,
                             aes(x = mod_var, y = coef,
                                 ymin = lwr, 
                                 ymax = upr)) +
  geom_hline(yintercept = 0,
             color = I(hsv(0/12, 7/12, 7/12)), 
             linetype = "dashed") + 
  geom_pointrange(size = 0.7) +
  facet_wrap(~ model, nrow = 1, scales = "free_x") +
  scale_y_continuous("Estimated DiD Effect on Dem. Vote Share",
                     limits = c(-0.03, 0.01)) +
  geom_text(data = txt.pos.nativism, 
            aes(label = var),
            size = axis.title.size - 12) +
  annotation_custom(b1) +
  annotate("text", x = 2, y = 0.01, 
           label = "2012 vs. 2016", 
           size = 4,
           hjust = 0.54) +
  annotation_custom(b2) +
  annotate("text", x = 5, y = 0.01, 
           label = "2012 vs. 2020", 
           size = 4,
           hjust = 0.46) +
  theme_bw() +  
  theme(legend.position = "none",
        axis.title.y = element_text(size = axis.title.size - 1,
                                    margin = margin(0, 20, 0, 0)),
        axis.text.y = element_text(size = axis.title.size - 1),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        plot.title = element_blank(),
        strip.text = element_text(size = axis.title.size - 1,
                                  face = "bold",
                                  margin = margin(0, 0, 5, 0)),
        strip.background = element_blank(),
        panel.grid = element_blank()
  ) 

# print
pdf(file = paste(MAIN_DIR, "Figure-3.pdf", sep = "/"), 
    width = 9, height = 4.5)
print(plot.coef.nativism)
dev.off()

setEPS()
cairo_ps(paste(MAIN_DIR, "Figure-3.eps", sep = "/"), 
         width = 9, height = 4.5)
print(plot.coef.nativism)
dev.off()


## Table D.5, Column 1 ---------------------------------------------------------
ddd.12.16.ho.bi.full <- felm(dem_share ~ 
                               anticorrupt + 
                               cn_ug_sqmi_cat_2012_bi_high:anticorrupt +
                               
                               anticorrupt:share_homeowner_county_2012_cat_medium + 
                               anticorrupt:share_homeowner_county_2012_cat_high + 
                               
                               anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium +
                               anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high +
                               
                               # controls
                               cn_g_sqmi_cat_2012_bi:anticorrupt +
                               ind_ug_sqmi_cat_2012_bi:anticorrupt +
                               
                               # demographics
                               share_pop_county_female * anticorrupt +
                               
                               share_pop_5_17_county * anticorrupt +
                               share_pop_18_24_county * anticorrupt + 
                               share_pop_25_34_county * anticorrupt +
                               share_pop_35_44_county * anticorrupt +
                               share_pop_45_54_county * anticorrupt +
                               share_pop_55_64_county * anticorrupt +
                               share_pop_65_county * anticorrupt + 
                               
                               share_pop_white_county * anticorrupt + 
                               share_pop_black_county * anticorrupt + 
                               share_pop_hispanic_county * anticorrupt + 
                               share_pop_asian_county * anticorrupt + 
                               share_pop_cn_county * anticorrupt +
                               
                               share_foreign_born_county * anticorrupt + 
                               
                               share_edu_ba_above_county * anticorrupt +
                               share_enroll_county * anticorrupt + 
                               ln_pop_density_county * anticorrupt + 
                               
                               effective_tax_county * anticorrupt + 
                               
                               ipw_county * anticorrupt + 
                               
                               employ_rate_county * anticorrupt + 
                               median_household_income_county * anticorrupt +
                               vacancy_county * anticorrupt 
                             | county_fips | 0 | county_fips, 
                             keepX = TRUE,
                             data = df.treat.1)
summary(ddd.12.16.ho.bi.full)
n_distinct(ddd.12.16.ho.bi.full$fe$county_fips)


## Table D.5, Column 2 ---------------------------------------------------------
ddd.12.20.ho.bi.full <- felm(dem_share ~ 
                               anticorrupt + 
                               cn_ug_sqmi_cat_2012_bi_high:anticorrupt +
                               
                               anticorrupt:share_homeowner_county_2012_cat_medium + 
                               anticorrupt:share_homeowner_county_2012_cat_high + 
                               
                               anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium +
                               anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high +
                               
                               # controls
                               cn_g_sqmi_cat_2012_bi:anticorrupt +
                               ind_ug_sqmi_cat_2012_bi:anticorrupt +
                               
                               # demographics
                               share_pop_county_female * anticorrupt +
                               
                               share_pop_5_17_county * anticorrupt +
                               share_pop_18_24_county * anticorrupt + 
                               share_pop_25_34_county * anticorrupt +
                               share_pop_35_44_county * anticorrupt +
                               share_pop_45_54_county * anticorrupt +
                               share_pop_55_64_county * anticorrupt +
                               share_pop_65_county * anticorrupt + 
                               
                               share_pop_white_county * anticorrupt + 
                               share_pop_black_county * anticorrupt + 
                               share_pop_hispanic_county * anticorrupt + 
                               share_pop_asian_county * anticorrupt + 
                               share_pop_cn_county * anticorrupt +
                               
                               share_foreign_born_county * anticorrupt + 
                               
                               share_edu_ba_above_county * anticorrupt +
                               share_enroll_county * anticorrupt + 
                               ln_pop_density_county * anticorrupt + 
                               
                               effective_tax_county * anticorrupt + 
                               
                               ipw_county * anticorrupt + 
                               
                               employ_rate_county * anticorrupt + 
                               median_household_income_county * anticorrupt + 
                               vacancy_county * anticorrupt 
                             | county_fips | 0 | county_fips, 
                             keepX = TRUE,
                             data = df.treat.2)
summary(ddd.12.20.ho.bi.full)
n_distinct(ddd.12.20.ho.bi.full$fe$county_fips)


## Table D.5, Column 3 ---------------------------------------------------------
ddd.12.16.v.bi.full <- felm(dem_share ~ 
                              anticorrupt + 
                              cn_ug_sqmi_cat_2012_bi_high:anticorrupt +
                              
                              anticorrupt:vacancy_county_2012_tri_medium + 
                              anticorrupt:vacancy_county_2012_tri_high + 
                              
                              anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium +
                              anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high +
                              
                              # demographics
                              share_pop_county_female * anticorrupt +
                              
                              share_pop_5_17_county * anticorrupt +
                              share_pop_18_24_county * anticorrupt + 
                              share_pop_25_34_county * anticorrupt +
                              share_pop_35_44_county * anticorrupt +
                              share_pop_45_54_county * anticorrupt +
                              share_pop_55_64_county * anticorrupt +
                              share_pop_65_county * anticorrupt + 
                              
                              share_pop_white_county * anticorrupt + 
                              share_pop_black_county * anticorrupt + 
                              share_pop_hispanic_county * anticorrupt + 
                              share_pop_asian_county * anticorrupt + 
                              share_pop_cn_county * anticorrupt + 
                              
                              share_foreign_born_county * anticorrupt + 
                              
                              share_edu_ba_above_county * anticorrupt +
                              share_enroll_county * anticorrupt + 
                              ln_pop_density_county * anticorrupt + 
                              
                              effective_tax_county * anticorrupt + 
                              
                              ipw_county * anticorrupt + 
                              
                              employ_rate_county * anticorrupt +
                              median_household_income_county * anticorrupt + 
                              vacancy_county * anticorrupt 
                            | county_fips | 0 | county_fips, 
                            keepX = TRUE,
                            data = df.treat.1)
summary(ddd.12.16.v.bi.full)
n_distinct(ddd.12.16.v.bi.full$fe$county_fips)


## Table D.5, Column 4 ---------------------------------------------------------
ddd.12.20.v.bi.full <- felm(dem_share ~ 
                              anticorrupt + 
                              cn_ug_sqmi_cat_2012_bi_high:anticorrupt +
                              
                              anticorrupt:vacancy_county_2012_tri_medium + 
                              anticorrupt:vacancy_county_2012_tri_high + 
                              
                              anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium +
                              anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high +
                              
                              # demographics
                              share_pop_county_female * anticorrupt +
                              
                              share_pop_5_17_county * anticorrupt +
                              share_pop_18_24_county * anticorrupt + 
                              share_pop_25_34_county * anticorrupt +
                              share_pop_35_44_county * anticorrupt +
                              share_pop_45_54_county * anticorrupt +
                              share_pop_55_64_county * anticorrupt +
                              share_pop_65_county * anticorrupt + 
                              
                              share_pop_white_county * anticorrupt + 
                              share_pop_black_county * anticorrupt + 
                              share_pop_hispanic_county * anticorrupt + 
                              share_pop_asian_county * anticorrupt + 
                              share_pop_cn_county * anticorrupt + 
                              
                              share_foreign_born_county * anticorrupt + 
                              
                              share_edu_ba_above_county * anticorrupt +
                              share_enroll_county * anticorrupt + 
                              ln_pop_density_county * anticorrupt + 
                              
                              effective_tax_county * anticorrupt + 
                              
                              ipw_county * anticorrupt + 
                              
                              employ_rate_county * anticorrupt + 
                              median_household_income_county * anticorrupt + 
                              vacancy_county * anticorrupt 
                            | county_fips | 0 | county_fips, 
                            keepX = TRUE,
                            data = df.treat.2)
summary(ddd.12.20.v.bi.full)
n_distinct(ddd.12.20.v.bi.full$fe$county_fips)


## Table D.5: Heterogeneous Treatment Effects: Pocketbook Voting Hypotheses ----
# set list of model outputs
did.pocket <- list(ddd.12.16.ho.bi.full,
                   ddd.12.20.ho.bi.full,
                   ddd.12.16.v.bi.full,
                   ddd.12.20.v.bi.full)

# set list of clustered SEs
did.pocket.cse <- list(ddd.12.16.ho.bi.full$cse,
                       ddd.12.20.ho.bi.full$cse,
                       ddd.12.16.v.bi.full$cse,
                       ddd.12.20.v.bi.full$cse)

# set model names
mod.names.pocket <- c("12-16-full",
                      "12-20-full",
                      "12-16-full",
                      "12-20-full")

# set variable order
var.order.pocket <- c("^anticorrupt$",
                      "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high$",
                      
                      "^anticorrupt:share_homeowner_county_2012_cat_medium$",
                      "^anticorrupt:share_homeowner_county_2012_cat_high$",
                      "^anticorrupt:vacancy_county_2012_tri_medium$",
                      "^anticorrupt:vacancy_county_2012_tri_high$",
                      
                      "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium$",
                      "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high$",
                      "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium$",
                      "^anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high$",
                      
                      # controls
                      "^anticorrupt:cn_g_sqmi_cat_2012_biCN-G-high$",
                      "^anticorrupt:ind_ug_sqmi_cat_2012_biIND-UG-high$",
                      
                      "^share_pop_county_female$",
                      "^share_pop_5_17_county$",
                      "^share_pop_18_24_county$",
                      "^share_pop_25_34_county$",
                      "^share_pop_35_44_county$",
                      "^share_pop_45_54_county$",
                      "^share_pop_55_64_county$",
                      "^share_pop_65_county$",
                      
                      "^share_pop_white_county$",
                      "^share_pop_black_county$",
                      "^share_pop_hispanic_county$",
                      "^share_pop_asian_county$",
                      "^share_pop_cn_county$",
                      
                      "^share_foreign_born_county$",
                      "^share_edu_ba_above_county$",
                      "^share_enroll_county$",
                      "^ln_pop_density_county$",
                      "^effective_tax_county$", 
                      "^ipw_county$",
                      "^employ_rate_county$", 
                      "^median_household_income_county$",
                      "^vacancy_county$",
                      
                      # trending effects of controls
                      "^anticorrupt:share_pop_county_female$",
                      "^anticorrupt:share_pop_5_17_county$",
                      "^anticorrupt:share_pop_18_24_county$",
                      "^anticorrupt:share_pop_25_34_county$",
                      "^anticorrupt:share_pop_35_44_county$",
                      "^anticorrupt:share_pop_45_54_county$",
                      "^anticorrupt:share_pop_55_64_county$",
                      "^anticorrupt:share_pop_65_county$",
                      
                      "^anticorrupt:share_pop_white_county$",
                      "^anticorrupt:share_pop_black_county$",
                      "^anticorrupt:share_pop_hispanic_county$",
                      "^anticorrupt:share_pop_asian_county$",
                      "^anticorrupt:share_pop_cn_county$",
                      
                      "^anticorrupt:share_foreign_born_county$",
                      "^anticorrupt:share_edu_ba_above_county$",
                      "^anticorrupt:share_enroll_county$",
                      "^anticorrupt:ln_pop_density_county$",
                      "^anticorrupt:effective_tax_county$", 
                      "^anticorrupt:ipw_county$",
                      "^anticorrupt:employ_rate_county$", 
                      "^anticorrupt:median_household_income_county$",
                      "^anticorrupt:vacancy_county$")


# set var label
var.label.pocket <- str_replace_all(var.order.pocket, "\\^", "")
var.label.pocket <- str_replace_all(var.label.pocket, "\\$", "")

# save
sink(file.path(MAIN_DIR, "Table-D5.txt"))
#sink(file.path(MAIN_DIR, "Table-D5.tex"))
stargazer(did.pocket,
          #type = "latex",
          type = "text",
          title = "Heterogeneous Treatment Effects: Pocketbook Voting Hypotheses",
          dep.var.labels = "Two-Party Democratic Vote Share",
          dep.var.labels.include = TRUE,
          column.labels = mod.names.pocket,
          order = var.order.pocket,
          covariate.labels = replaceVarName(var.vec = var.label.pocket,
                                            var.df = var.df),
          add.lines = list(c("Fixed Effects: County",
                             rep("Yes", 4)),
                           c("Number of Counties",
                             formatC(n_distinct(ddd.12.16.ho.bi.full$fe$county_fips), big.mark = ","),
                             formatC(n_distinct(ddd.12.20.ho.bi.full$fe$county_fips), big.mark = ","),
                             formatC(n_distinct(ddd.12.16.v.bi.full$fe$county_fips), big.mark = ","),
                             formatC(n_distinct(ddd.12.20.v.bi.full$fe$county_fips), big.mark = ",")
                           )),
          label = "tb:did-pres-moderation-pocket",
          single.row = FALSE,
          se = did.pocket.cse,
          df = FALSE,
          font.size = "tiny",
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          star.char = c("+", "*", "**", "***"),
          notes.append = FALSE)
sink()


## Figure 4: Effect Heterogeneity: Pocketbook Voting Mechanism -----------------
# extract DiD results in list
fm.list.pocket <- list("Initial Homeownership Rate\n" = summary(glht(ddd.12.16.ho.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                       "Initial Homeownership Rate\n" = summary(glht(ddd.12.16.ho.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium = 0"))),
                       "Initial Homeownership Rate\n" = summary(glht(ddd.12.16.ho.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high = 0"))),
                       
                       "Initial Homeownership Rate\n" = summary(glht(ddd.12.20.ho.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                       "Initial Homeownership Rate\n" = summary(glht(ddd.12.20.ho.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium = 0"))),
                       "Initial Homeownership Rate\n" = summary(glht(ddd.12.20.ho.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high = 0"))),
                       
                       "Initial Vacancy Rate\n" = summary(glht(ddd.12.16.v.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                       "Initial Vacancy Rate\n" = summary(glht(ddd.12.16.v.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium = 0"))),
                       "Initial Vacancy Rate\n" = summary(glht(ddd.12.16.v.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high = 0"))),
                       
                       "Initial Vacancy Rate\n" = summary(glht(ddd.12.20.v.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                       "Initial Vacancy Rate\n" = summary(glht(ddd.12.20.v.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium = 0"))),
                       "Initial Vacancy Rate\n" = summary(glht(ddd.12.20.v.bi.full, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high = 0")))
)

# extract key coefs and clustered SEs
coef.df.pocket <- map_df(fm.list.pocket, function(x){ 
  out <- tibble(var = names(coef(x)),
                coef = coef(x),
                cse = sqrt(diag(vcov(x))))
}, .id = "model") 

# add variables
var.name.pocket <- c("Low",
                     "Medium",
                     "High",
                     "Low",
                     "Medium",
                     "High")

var.period.pocket <- c(rep("2012 vs. 2016", 3),
                       rep("2012 vs. 2020", 3))

# combine
coef.df.pocket <- cbind(coef.df.pocket, 
                        var.name.pocket,
                        var.period.pocket)

# set alpha level
alpha <- 0.05

# compute C.I.
coef.df.pocket <- coef.df.pocket %>%
  mutate(model = factor(model, 
                        levels = c("Initial Homeownership Rate\n",
                                   "Initial Vacancy Rate\n")),
         multiplier = qnorm(1 - alpha / 2),
         coef = as.numeric(coef),
         cse = as.numeric(cse),
         lwr = coef - multiplier * cse,
         upr = coef + multiplier * cse,
         mod_var = paste(model, var.name.pocket, var.period.pocket, sep = "-"))

# set variable order
coef.df.pocket <- coef.df.pocket %>%
  mutate(mod_var = factor(mod_var, 
                          levels = coef.df.pocket$mod_var))

# set text df for annotations
txt.shift.pocket <- 0.002
txt.pos.pocket <- tibble(model = coef.df.pocket$model,
                         var = as.character(coef.df.pocket$var.name),
                         mod_var = coef.df.pocket$mod_var,
                         coef = coef.df.pocket$lwr - txt.shift.pocket,
                         lwr = coef.df.pocket$lwr,
                         upr = coef.df.pocket$upr)

# set parameters
axis.title.size <- 16

# add curly brackets
bracketsGrob <- function(...){
  l <- list(...)
  e <- new.env()
  e$l <- l
  grid:::recordGrob(  {
    do.call(grid.brackets, l)
  }, e)
}

brack.y <- 0.87
b1 <- bracketsGrob(0.05, brack.y, 0.45, brack.y, 
                   h = 0.05, lwd = 1,
                   col = "black")
b2 <- bracketsGrob(0.55, brack.y, 0.95, brack.y, 
                   h = 0.05, lwd = 1,
                   col = "black")

# plot
plot.coef.pocket <- ggplot(data = coef.df.pocket,
                           aes(x = mod_var, y = coef,
                               ymin = lwr, 
                               ymax = upr)) +
  geom_hline(yintercept = 0,
             color = I(hsv(0/12, 7/12, 7/12)), 
             linetype = "dashed") + 
  geom_pointrange(size = 0.7) +
  facet_wrap(~ model, nrow = 1, scales = "free_x") +
  scale_y_continuous("Estimated DiD Effect on Dem. Vote Share",
                     limits = c(-0.03, 0.01)) +
  geom_text(data = txt.pos.pocket, 
            aes(label = var),
            size = axis.title.size - 12) +
  annotation_custom(b1) +
  annotate("text", x = 2, y = 0.01, 
           label = "2012 vs. 2016", 
           size = 4,
           hjust = 0.54) +
  annotation_custom(b2) +
  annotate("text", x = 5, y = 0.01, 
           label = "2012 vs. 2020", 
           size = 4,
           hjust = 0.46) +
  theme_bw() +  
  theme(legend.position = "none",
        axis.title.y = element_text(size = axis.title.size - 1,
                                    margin = margin(0, 20, 0, 0)),
        axis.text.y = element_text(size = axis.title.size - 1),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        plot.title = element_blank(),
        strip.text = element_text(size = axis.title.size - 1,
                                  face = "bold",
                                  margin = margin(0, 0, 5, 0)),
        strip.background = element_blank(),
        panel.grid = element_blank()
  ) 

# print
pdf(file = paste(MAIN_DIR, "Figure-4.pdf", sep = "/"), 
    width = 9, height = 4.5)
print(plot.coef.pocket)
dev.off()


setEPS()
cairo_ps(paste(MAIN_DIR, "Figure-4.eps", sep = "/"), 
         width = 9, height = 4.5)
print(plot.coef.pocket)
dev.off()


## Figure D.2: Heterogeneous Treatment Effects: Joint Model Results ------------
# 2012 vs. 2016
ddd.12.16.all <- felm(dem_share ~ 
                        anticorrupt + 
                        
                        # CN FREI exposure
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high + 
                        
                        # white share
                        anticorrupt:share_white_2012_cat_medium + 
                        anticorrupt:share_white_2012_cat_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high +
                        
                        # BA above share
                        anticorrupt:share_edu_ba_above_county_2012_tri_medium + 
                        anticorrupt:share_edu_ba_above_county_2012_tri_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high +
                        
                        # homeownership rate
                        anticorrupt:share_homeowner_county_2012_cat_medium + 
                        anticorrupt:share_homeowner_county_2012_cat_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high +
                        
                        # vacancy rate
                        anticorrupt:vacancy_county_2012_tri_medium + 
                        anticorrupt:vacancy_county_2012_tri_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high +
                        
                        # controls
                        cn_g_sqmi_cat_2012_bi:anticorrupt +
                        ind_ug_sqmi_cat_2012_bi:anticorrupt +
                        
                        # demographics
                        share_pop_county_female * anticorrupt +
                        
                        share_pop_5_17_county * anticorrupt +
                        share_pop_18_24_county * anticorrupt + 
                        share_pop_25_34_county * anticorrupt +
                        share_pop_35_44_county * anticorrupt +
                        share_pop_45_54_county * anticorrupt +
                        share_pop_55_64_county * anticorrupt +
                        share_pop_65_county * anticorrupt + 
                        
                        share_pop_white_county * anticorrupt + 
                        share_pop_black_county * anticorrupt + 
                        share_pop_hispanic_county * anticorrupt + 
                        share_pop_asian_county * anticorrupt + 
                        share_pop_cn_county * anticorrupt + 
                        
                        share_foreign_born_county * anticorrupt + 
                        
                        share_edu_ba_above_county * anticorrupt +
                        share_enroll_county * anticorrupt + 
                        ln_pop_density_county * anticorrupt + 
                        
                        effective_tax_county * anticorrupt + 
                        
                        ipw_county * anticorrupt + 
                        
                        employ_rate_county * anticorrupt + 
                        median_household_income_county * anticorrupt + 
                        vacancy_county * anticorrupt 
                      | county_fips | 0 | county_fips, 
                      keepX = TRUE,
                      data = df.treat.1)
summary(ddd.12.16.all)


# 2012 vs. 2020
ddd.12.20.all <- felm(dem_share ~ 
                        anticorrupt + 
                        
                        # CN FREI exposure
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high + 
                        
                        # white share
                        anticorrupt:share_white_2012_cat_medium + 
                        anticorrupt:share_white_2012_cat_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high +
                        
                        # BA above share
                        anticorrupt:share_edu_ba_above_county_2012_tri_medium + 
                        anticorrupt:share_edu_ba_above_county_2012_tri_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high +
                        
                        # homeownership rate
                        anticorrupt:share_homeowner_county_2012_cat_medium + 
                        anticorrupt:share_homeowner_county_2012_cat_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high +
                        
                        # vacancy rate
                        anticorrupt:vacancy_county_2012_tri_medium + 
                        anticorrupt:vacancy_county_2012_tri_high + 
                        
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium +
                        anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high +
                        
                        # controls
                        cn_g_sqmi_cat_2012_bi:anticorrupt +
                        ind_ug_sqmi_cat_2012_bi:anticorrupt +
                        
                        # demographics
                        share_pop_county_female * anticorrupt +
                        
                        share_pop_5_17_county * anticorrupt +
                        share_pop_18_24_county * anticorrupt + 
                        share_pop_25_34_county * anticorrupt +
                        share_pop_35_44_county * anticorrupt +
                        share_pop_45_54_county * anticorrupt +
                        share_pop_55_64_county * anticorrupt +
                        share_pop_65_county * anticorrupt + 
                        
                        share_pop_white_county * anticorrupt + 
                        share_pop_black_county * anticorrupt + 
                        share_pop_hispanic_county * anticorrupt + 
                        share_pop_asian_county * anticorrupt + 
                        share_pop_cn_county * anticorrupt + 
                        
                        share_foreign_born_county * anticorrupt + 
                        
                        share_edu_ba_above_county * anticorrupt +
                        share_enroll_county * anticorrupt + 
                        ln_pop_density_county * anticorrupt + 
                        
                        effective_tax_county * anticorrupt + 
                        
                        ipw_county * anticorrupt + 
                        
                        employ_rate_county * anticorrupt + 
                        median_household_income_county * anticorrupt + 
                        vacancy_county * anticorrupt 
                      | county_fips | 0 | county_fips, 
                      keepX = TRUE,
                      data = df.treat.2)
summary(ddd.12.20.all)


# extract DiD results in list
fm.list.combined <- list("Initial White Population\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial White Population\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium = 0"))),
                         "Initial White Population\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high = 0"))),
                         
                         "Initial White Population\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial White Population\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_medium = 0"))),
                         "Initial White Population\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_white_2012_cat_high = 0"))),
                         
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium = 0"))),
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high = 0"))),
                         
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_medium = 0"))),
                         "Initial College-Educated Population\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_edu_ba_above_county_2012_tri_high = 0"))),
                         "Initial Homeownership Rate\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial Homeownership Rate\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium = 0"))),
                         "Initial Homeownership Rate\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high = 0"))),
                         
                         "Initial Homeownership Rate\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial Homeownership Rate\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_medium = 0"))),
                         "Initial Homeownership Rate\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:share_homeowner_county_2012_cat_high = 0"))),
                         
                         "Initial Vacancy Rate\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial Vacancy Rate\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium = 0"))),
                         "Initial Vacancy Rate\n" = summary(glht(ddd.12.16.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high = 0"))),
                         
                         "Initial Vacancy Rate\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high = 0"))),
                         "Initial Vacancy Rate\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_medium = 0"))),
                         "Initial Vacancy Rate\n" = summary(glht(ddd.12.20.all, linfct = c("anticorrupt:cn_ug_sqmi_cat_2012_bi_high + anticorrupt:cn_ug_sqmi_cat_2012_bi_high:vacancy_county_2012_tri_high = 0")))
)

# extract key coefs and clustered SEs
coef.df.combined <- map_df(fm.list.combined, function(x){ 
  out <- tibble(var = names(coef(x)),
                coef = coef(x),
                cse = sqrt(diag(vcov(x))))
}, .id = "model") 

# add variables
var.name.combined <- c("Small",
                       "Medium",
                       "Large",
                       "Small",
                       "Medium",
                       "Large")

var.period.combined <- c(rep("2012 vs. 2016", 3),
                         rep("2012 vs. 2020", 3),
                         rep("2012 vs. 2016", 3),
                         rep("2012 vs. 2020", 3),
                         rep("2012 vs. 2016", 3),
                         rep("2012 vs. 2020", 3),
                         rep("2012 vs. 2016", 3),
                         rep("2012 vs. 2020", 3))

# combine
coef.df.combined <- cbind(coef.df.combined, 
                          var.name.combined,
                          var.period.combined)

# set alpha level
alpha <- 0.05

# compute C.I.
coef.df.combined <- coef.df.combined %>%
  mutate(model = factor(model, 
                        levels = c("Initial White Population\n",
                                   "Initial College-Educated Population\n",
                                   "Initial Homeownership Rate\n",
                                   "Initial Vacancy Rate\n"
                        )),
         multiplier = qnorm(1 - alpha / 2),
         coef = as.numeric(coef),
         cse = as.numeric(cse),
         lwr = coef - multiplier * cse,
         upr = coef + multiplier * cse,
         mod_var = paste(model, var.name.combined, var.period.combined, sep = "-"))

# set variable order
coef.df.combined <- coef.df.combined %>%
  mutate(mod_var = factor(mod_var, 
                          levels = coef.df.combined$mod_var))

# set text df for annotations
txt.shift.combined <- 0.002
txt.pos.combined <- tibble(model = coef.df.combined$model,
                           var = as.character(coef.df.combined$var.name),
                           mod_var = coef.df.combined$mod_var,
                           coef = coef.df.combined$lwr - txt.shift.combined,
                           lwr = coef.df.combined$lwr,
                           upr = coef.df.combined$upr)

# set parameters
axis.title.size <- 16

# # add curly brackets
bracketsGrob <- function(...){
  l <- list(...)
  e <- new.env()
  e$l <- l
  grid:::recordGrob(  {
    do.call(grid.brackets, l)
  }, e)
}

brack.y <- 0.87
b1 <- bracketsGrob(0.05, brack.y, 0.45, brack.y, 
                   h = 0.05, lwd = 1,
                   col = "black")
b2 <- bracketsGrob(0.55, brack.y, 0.95, brack.y, 
                   h = 0.05, lwd = 1,
                   col = "black")

# plot
plot.coef.combined <- ggplot(data = coef.df.combined,
                             aes(x = mod_var, y = coef,
                                 ymin = lwr, 
                                 ymax = upr)) +
  geom_hline(yintercept = 0,
             color = I(hsv(0/12, 7/12, 7/12)), 
             linetype = "dashed") + 
  geom_pointrange(size = 0.7) +
  facet_wrap(~ model, nrow = 1, scales = "free_x") +
  scale_y_continuous("Estimated DiD Effect on Dem. Vote Share",
                     limits = c(-0.04, 0.03)) +
  geom_text(data = txt.pos.combined, 
            aes(label = var),
            size = axis.title.size - 12) +
  annotation_custom(b1) +
  annotate("text", x = 2, y = 0.03, 
           label = "2012 vs. 2016", 
           size = 4,
           hjust = 0.54) +
  annotation_custom(b2) +
  annotate("text", x = 5, y = 0.03, 
           label = "2012 vs. 2020", 
           size = 4,
           hjust = 0.46) +
  theme_bw() +  
  theme(legend.position = "none",
        axis.title.y = element_text(size = axis.title.size - 1,
                                    margin = margin(0, 20, 0, 0)),
        axis.text.y = element_text(size = axis.title.size - 1),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        plot.title = element_blank(),
        strip.text = element_text(size = axis.title.size - 1,
                                  face = "bold",
                                  margin = margin(0, 0, 5, 0)),
        strip.background = element_blank(),
        panel.grid = element_blank()) 


# print
pdf(file = paste(MAIN_DIR, "Figure-D2.pdf", sep = "/"), 
    width = 18, height = 4.5)
print(plot.coef.combined)
dev.off()

